#8/12/2015
#replication materials for Spirling "Democratization and Linguistic Complexity", JOP

#Figure 6 is estd bhat on the session-by-session regressions.
# This script also does the structural break test mentioned in text

rm(list=ls())

#load the data
#setwd("C:/Users/as9934/Dropbox/complexity/August2015/JOP_replication_data/")
load("bigframe.rdata")



#################################
#### Figure 6 ###################
#################################

sessions <- as.character( unique(big.frame$year.dummy) )

#matrix to store results
coef.frame <- data.frame(sess=unique(big.frame$year.dummy), est = as.numeric(NA), se = as.numeric(NA),
                         lo=as.numeric(NA), hi=as.numeric(NA), pval=as.numeric(NA))


#clus SEs function (courtesy of Drew Dimmery)
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }



for(i in 1:length(sessions)){
  
  
  dat <- big.frame[big.frame$year.dummy==sessions[i],]
  dat <- dat[complete.cases(dat),]
  
  
  #main regression
  mod <- lm(FK_read.ease~ party + cabinet+ word.count + competitiveness, data= dat)
  
  
  
  #fix SEs
  # apply the 'cl' function , cluster on mp
  mod.clus.SE <-cl(dat, mod, dat$mp_code)
  
  
  coefs <- mod.clus.SE
  coef.frame[i,c(2,3)] <- as.numeric(coefs["cabinet1",c(1,2)])
  
  coef.frame[i,4] <- coefs["cabinet1",c(1)] + ( 1.96*as.numeric(coefs["cabinet1",2]))
  coef.frame[i,5] <- coefs["cabinet1",c(1)] -( 1.96*as.numeric(coefs["cabinet1",2]))
  coef.frame[i,6] <- coefs["cabinet1",4]
  
  cat("\n done",sessions[i],"\n")
  
}

ylimup <-  max(c(coef.frame$lo, coef.frame$hi))
ylimdown <- min(c(coef.frame$lo, coef.frame$hi))




par(bg='cornsilk1')
plot(1:nrow(coef.frame), coef.frame$est, type="p", pch=16, ylim=c(ylimdown, ylimup), axes=F, xlab="",
     ylab="")


mtext(expression(hat(beta)), side=2, line=2, cex=1.5, srt=90 )

par(las=1)
axis(2)
axis(1, at=1:nrow(coef.frame), labels=coef.frame$sess, cex.axis=.7, las=1)
box()

arrows(1:nrow(coef.frame), coef.frame$lo, 1:nrow(coef.frame), coef.frame$hi, code=3, angle=90, length=.05)
abline(h=0, col="red", lwd=2)  

#do structural break test
require(strucchange)
y <- coef.frame$est
x <- 1:nrow(coef.frame)
print( breakpoints(y~x) )